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This work explores the use of a tree tensor network ansatz to simulate the ground state of a 
local Hamiltonian on a two-dimensional lattice. By exploiting the entropic area law, the tree tensor 
network ansatz seems to produce quasi-exact results in systems with sizes well beyond the reach of 
exact diagonalisation techniques. We describe an algorithm to approximate the ground state of a 
local Hamiltonian on a L x L lattice with the topology of a torus. Accurate results are obtained 
for L = {4,6,8}, whereas approximate results are obtained for larger lattices. As an application 
of the approach, we analyse the scaling of the ground state entanglement entropy at the quantum 
critical point of the model. We confirm the presence of a positive additive constant to the area law 
for half a torus. We also find a logarithmic additive correction to the entropic area law for a square 
block. The single copy entanglement for half a torus reveals similar corrections to the area law with 
a further term proportional to 1/L. 
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I. INTRODUCTION 

The numerical study of many-body quantum systems 
is a challenging task. The exponential growth of the di- 
mension of the Hilbert space with the size of the system 
implies that exact diagonalisation techniques can only be 
applied to address small lattice systems [Ij |2l |3] . Quan- 
tum Monte Carlo sampling offers a valuable route to the 
study of larger lattices, although systems of frustrated 
quantum spins or interacting fermions cannot be anal- 
ysed due to the so called sign problem. 

In two spatial dimensions, the use of a tensor net- 
work ansatz, such as the tensor product state or pro- 
jected entangled pair state (PEPS) [llSlll|71IHlli[inj[II] 
and the multi-scale entanglement renormalisation ansatz 
(MERA) [121 El El E], has opened a very promising 
alternative path to investigating ground state properties 
of arbitrarily large lattice systems. The key of these ap- 
proaches is the abihty of the TPS, PEPS and MERA to 
reproduce the scaling of the entanglement in the ground 
state, as given by the entropic area law. 

In this work we explore the use of yet another tensor 
network variational ansatz, namely a tree tensor network 
(TTN) [16l[17lElEll20l|2ll|221|23], to simulate the 
ground state of local 2D lattice systems. This very sim- 
ple ansatz is inspired on the original real space Renor- 
malisation Group ideas of Kadanoff, Migdal and Wilson 

The present approach is both motivated and limited 
by the area law for the entanglement entropy [501 1311 1311 
1331 1311 13S]- The "area law" is a conjectured property of 
the entanglement entropy of certain ground states of lo- 
cal Hamiltonians. It asserts that the entropy of a region 
of the system is proportional to the size of its boundary 
rather than proportional to its volume. Direct calcula- 
tions of the entanglement entropy have provided evidence 
of the validity of the area law for a large class of sys- 



tems [3ni [311 132 133] ■ A complete characterization of the 
Hamiltonians whose ground state obeys the area law is 
still missing. 

On the one hand, by exploiting the area law a TTN 
can be used to address small 2D lattices with sizes well 
beyond the reach of exact diagonalisation techniques. 
Specifically, the cost of simulating a lattice of L x L sites 
grows as exp(L) instead of exp(L^). Thus, the TTN ap- 
proach is useful to investigate small 2D quantum sys- 
tems and to study larger systems with finite size scaling 
techniques. It is also particularly suitable to investigate 
ground state entropies. 

On the other hand, the exp(L) cost due to the en- 
tropic area law still sets a severe limit to the system sizes 
a TTN can describe and the present approach simply 
cannot compete with the PEPS and MERA algorithms 
[1 [SI Uni HH UD for large systems. However, the TTN is 
also of interest in the context of developing these more 
advanced, scalable algorithms. This is due both to its 
simplicity and to its direct connection to ground state en- 
tanglement properties, on which the scalable algorithms 
are also based. As a matter of fact, the TTN approach 
described in this work was initially developed as an aux- 
iliary tool to help in the design of the MERA [T^ . 

The present approach bears important similarities with 
White's density matrix renormalisation group (DMRG) 
[3S1I37113H] (for a review see e.g. P5^) when applied to 2D 
lattice systems [IHlllIljiailS! 44, 45, 46 . Roughly speak- 
ing, it can be regarded as a DMRG approach where the 
matrix product state has been replaced with a TTN. This 
replacement has both advantages and disadvantages. Its 
weakest point is an increase in computational cost. How- 
ever, a TTN greatly improves the connectivity between 
lattice sites, possibly resulting in faster convergence and 
better correlation functions (e.g. on a torus). Extract- 
ing certain entropies from the TTN, say the entropy of 
one quarter of the lattice, is straightforward. Finally, the 
algorithm can be very simply implemented. 



The results are organised in several sections. In Sect. 
In] we describe the TTN for 2D lattices and motivate its 
use in terms of the entropic area law. In Sect. |III| we 
explain how to compute the expectation value of local 
operators, two-point correlation functions, fidelities and 
block entropies. Then in Sect. |IV| we describe an algo- 
rithm to approximate ground states with a TTN. This 
algorithm is tested in Sect. |V]by addressing the quan- 
tum Ising model with transverse magnetic field on a torus 
made of L x L sites. Quasi exact results are obtained for 
lattices of linear size L = {4, 6, 8}, whereas approximate 
results are obtained for L — {10,16,32}, we also check 
the exponential cost of a TTN representation and intro- 
duce a possible estimate of the error induced by the finite 
amount of computational resources. 

In Sect. VI we turn our attention to the computation 
of ground state entanglement at the critical point. We 
present results for L = {4, 6, 8, 10} for both entangle- 
ment entropy and single copy entanglement. This allows 
us to investigate the form of their finite size scaling. Sev- 
eral authors |47l|48l|49l[50l|5Tl[52l|53IMl[55l|5^ have 
predicted the presence of corrections to the area law in- 
cluding a universal and positive constant term for half a 
torus |Sni EU m]. Our results reproduce this term, as 
well as a logarithmic correction for a quarter of a torus 
1521 ESI [SB] . Numerical estimates of the coefficients 
for all the corrections are presented in Eqs. |4r] |42] |45] 
I49II5T] We conclude with a discussion of the results in 
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FIG. 1: Example of TTN for a 2 x 2 lattice and a 4 x 4 lattice. 
Notice (right) that the TTN for a 2D lattice can always be 
represented as a planar graph, with the leaves or physical 
indices ordered on a line. The tensors labelled with Wi are 
isometric tensors. They act locally by projecting the ground 
state onto its local support with dimension Xi (see section JTl] 
for further explanation). 



II. TREE TENSOR NETWORK ANSATZ 

In this section we introduce the variational ansatz used 
throughout the manuscript and justify its applicability in 
terms of the area law for entanglement entropy. 

A. Isometric tree tensor network 

We consider a square lattice £ made of N — Lx L sites, 
where each site is described by a local Hilbert space V 
of finite dimension d. Our goal is to represent a pure 
state l^-) e V^^ of the lattice C. Most of the time, \^) 
will correspond to the ground state I^'gs) of some local 
Hamiltonian H defined on C. 

A generic state \'^) g can always be expanded as 

d d d 

ii=l 12=1 = ! 

where the coefficients Ti^i^...i^ are complex numbers 
and the vectors {|ls), |2s), • • • , Ms)} denote a local basis 
on site s € £. We refer to the index is that labels a local 
basis for site s (z^ = 1 , • • • ,d) as a physical index. 

In this work we further expand the tensor of coeffi- 
cients Ti^i^...i^ in Eq. [T] using a TTN. As shown in Fig. 
[ijfor lattices of 2 x 2 and 4x4 sites, a TTN decomposition 
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FIG. 2: (i) Diagrammatic representation of three types of 
isometric tensors in the TTN for a 4 x 4 lattice in Fig. [T| 
(ii) Graphical representation of the constraints in Eqs. |3j[5] 
fulfilled by the isometric tensors. 



consists of a collection of tensors w that have both bond 
indices and physical indices. The tensors are intercon- 
nected by the bond indices according to a tree pattern. 
The N physical indices correspond to the leaves of the 
tree. Upon summing over all the bond indices, the TTN 
produces the complex coefficients Ti^i^...i^ of Eq. [l] 

The tensors in the TTN will be constrained to be iso- 
metric, in the following sense. As shown in Fig. |2] for 
the 4x4 case of Fig. [T] each tensor w in a TTN has at 
most one upper leg/index a and some number p of lower 
indices/legs • • • , /3p, so that its entries read {w)% ..o . 
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Then we impose that 

E HX---pMfJ'''''' = ^^^'- (2) 

/3i--/3p 

For the sake of clarity, throughout the paper we use di- 
agrams to represent tensor networks as well as tensor 
manipulations. For instance, the constraints for the ten- 
sors Wi, W2 and of the TTN of Fig. [T] for a 4 x 4 
lattice, namely 

E MlpMi^lt''^''''' = (3) 

/3i/32/33/34 

E(«^2)^,ft(u;J)f;''^ = (4) 

Pi 132 

T^M^Mf'^' = 1: (5) 

are represented as diagrams in Fig. |2jii). We refer to a 
tensor w that fulfils Eq. [2] as an isometry. As we will see 
in Sects. |III| and [TVl the use of isometries simplifies the 
manipulations necessary to compute sxpectation values 
of local operators and the spectrum of reduced density 
matrices, as well as to optimise the TTN. The isomet- 
ric character of the tensors can also be seen to prevent 
numerical instability during the simulations. 



B. Coarse-graining of the lattice 

An intuitive interpretation of the use of a TTN to rep- 
resent a state \'^) can be obtained in terms of a coarse- 
graining transformation for the lattice C Notice that 
the isometries w in Fig. [l] are organised in layers. The 
bond indices between two layers can be interpreted as 
defining the sites of an effective lattice. In other words, 
the TTN defines a sequence of increasingly coarser lat- 
tices {£oj ^1, ■ ■ ■ , ^T-i}, where Co = C and each site of 
lattice Cr is defined in terms of several sites of Cr-i by 
means of an isometry Wr, see Fig. [3] 

In this picture, a site of the lattice £r effectively corre- 
sponds to some number rir of sites of the original lattice 
Cq. For instance, each of the two sites of £2 in Fig- [3] 
corresponds to 8 sites of Cq. Similarly, each site of lattice 
£1 corresponds to 4 sites of Lq- 



C. Entropic area law 

In using a TTN to represent a generic state \'^), the 
top tensor wt must already depend on c?^ coefficients. It 
is then unclear that the use of a TTN has any computa- 
tional advantage with respect to directly dealing with all 
the coefficients Ti^i^...ij^ in Eq. [l] However, ground 
states I^'gs) of local Hamiltonians are known to often ex- 
hibit a so-called entropic area law [501 1311 1321 ESI IMl 131] 




FIG. 3; The isometric TTN of Fig. [T] for a 4 x 4 lattice 
Co is associated with a coarse-graining transformation that 
generates a sequence of increasingly coarse-grained lattices 
£1, £2 and £3. Notice that in this example we have added 
an extra index to the top isometry wz, corresponding to the 
single site of an extra top lattice £3, which we can use to 
encode in the TTN a whole subspace of V®^ instead of a 
single state j^*). 



and this property might lead to a reduction in compu- 
tational costs when expressing the state (or an accurate 
approximation of it) as a TTN. 

Let us introduce the reduced density matrix p for a 
block A of contiguous sites of £ as 

/^trslWHE^'-l^-X*^^!' (6) 

where B are all the sites of £ outside the block A and 
Pa are the eigenvalues of p"^ (that is, p^Y^"^) — pa\^a))- 
Then the entropy S{A) of block A is defined as 

5(^) = -tr(p^l0gp^) = -EPal0gPa. (7) 

a 

This entropy measures the amount of entanglement be- 
tween the block A and the rest B of the lattice £, and it is 
also known as entanglement entropy 57J. For a generic 
state, the entropy of a block A is proportional to the 
number n{A) of sites in A (provided n{A) < N/2), that 
is 

S{A) w n{A) log d (generic) (8) 

For instance, the entropy of a block made of I x I sites is 
proportional to P and, correspondingly, the effective di- 
mension X required to describe the block ^ is exponential 



Throughout the paper we use a number of different subscripts 
to denote different effective dimensions x- For instance, Xixi: 
x(j4) and Xl/2 refer, respectively, to the effective dimension for 
a block of Z X ( sites, a block A and one half of the lattice. The 
specific meaning should be clear from the context. 



4 



in f 



Siy^i « l^logd, 



Xlxl 



(generic) 



(9) 



If instead the entropy of the block A grows propor- 
tional to the number of sites of the boundary of A, de- 
noted by we say that the state j^*) fulfils an en- 
tropic " area law" 



S'(A)«ciE(A) (area law) 



(10) 



where Ci is some constant. For instance, for the above 
block of I X I sites, the entropy is only proportional to 
I. Accordingly, the dimension x required to effectively 
describe the block may grow markedly less with I than 
in the generic case, 

Sixifnicil, X;xi > exp(4ci/), (area law) (11) 

where the lower bound for x is obtained by exponenti- 
ating the entropy and is saturated by a flat probability 
distribution pa ~ l/x, a = 1, • • • , X- 

Eq. [TT]is our main justification for attempting to de- 
scribe ground states of local 2D Hamiltonians using a 
TTN. It suggests that it might be possible to accurately 
approximate a ground state I^'gs) that fulfils the area law 
of Eq. [To] by using a number of coefficients that scales 
with the linear size L of the lattice £ only as 0(exp(L)), 
instead of 0(exp(i^)) as is the case for a generic state. 
In other words, ground states of local Hamiltonians are 
typically less entangled than generic states, and we might 
be able to exploit this fact computationally. 

A simple example of ground state that fuUfills the 
boundary law is a valence bond crystal (as e.g. it is 
the ground state of the AKLT model [5H]) made of short- 
range dimers. Each dimer crossing the boundaries of a re- 
gion A contributes a constant amount to the entropy and 
the number of such dimers is proportional to the size of 
the boundary of the region A. This roughly corresponds 
to the presence, in the spectrum of the reduced density 
matrix, of a number of significant eigenvalues that grows 
exponentially with the size of the boundary. 



D. Plane, cylinder and torus 

Let us now assume that the entropic boundary law in 
Eq. [To] translates into an effective site dimension given 

by 



x(A)«exp(5(A))«exp(cia(A)), 



(12) 



and let us explore the implications that this expression 
would have on the ability of a TTN to encode ground 
states. 

For this purpose, let us consider the (interacting) 
boundaries, denoted S1/2, and Si/g, of blocks that 
consists, respectively, of one half, one fourth and one 
eighth of a L X L lattice C These boundaries depend 
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TABLE I: Length of the boundaries of different subregions 
of a 2D lattice system for several choices of topology of the 
whole system. 
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FIG. 4: Interacting boundaries S1/2, and Ei/g corre- 

sponding to one half, one quarter and one eighth of a lattice 
for three different choices of boundary conditions. 



on the topology of the interactions of H on C, and for 
the plane, cylinder and torus are presented in table [l] (see 
also Fig. [4]). 

From table [I] and Eq. [12 one can obtain the dimension 
X of the sites of the most coarse-grained lattices Ct-i, 
Ct-2, J^T-Zi and the size of the isometrics at the upper 
layers of the TTN, which is what dominates the com- 
putational cost of the approach. The table shows that 
ground states on a torus are more entangled (e.g. the 
blocks have more interacting boundary, or entropy), and 
therefore computationally more demanding, than on a 
plane or cylinder. In this work we shall concentrate on 
the torus, with the understanding that a similar analysis 
can also be conducted for the other cases. [In particular, 
as it is easy to anticipate, given the same computational 
costs, larger systems can be addressed in the cases of 
plane and cylinder interaction topologies.] 



E. TTN ansatz on the torus 

From now on we consider a L x L lattice C on the 
torus. In this case, Xi/2 ~ Xi/4 ~ exp(ci2L) are the 
largest effective site dimensions. The top isometry wt 
depends on X1/2 ~ exp(4ciL) parameters, whereas each 
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isometry wt-i depends on X1/2X1/4 ~ exp(6ciL) param- 
eters. Isometries at lower layers of the TTN can be seen 
to depend on less parameters. 

Based on these observations, our TTN ansatz for the 
ground state of an L x L lattice with torus topology and 
site dimension d = 2 (e.g. spin-^ model) will invari- 
ably consists of a top isometry wt and two isometries 
wt-1 with bond dimension x on all their indices. Then, 
depending on the size L and other considerations, the 
TTN will be completed in two possible ways. For small 
L (L < 8 in the examples of section [V]), a single extra 
layer of isometries will be considered, where each isome- 
try maps N/ 4 sites of C directly into one site of Ct-2 ■ For 
larger lattices, it is computationally favourable to com- 
plete the TTN with at least two more layers of isometries, 
see Fig. [5j 

Because the isometries wt-i are, by far, the largest 
tensors in the TTN, the memory required to store the 
ansatz is a function of the size of wt^i, namely 



TTN for 6 X 6 Lattice 



Memory 



X 



(large x regime) (13) 



where, unless otherwise specified, from now on x refers 
to the effective dimension used for both one half and one 
quarter of the lattice, x = X1/2 — Xi/4- 



F. Nested Schmidt decompositions 

It is instructive to relate the TTN ansatz with the 
Schmidt decomposition of the state \'^) it represents. 

Recall that given a bipartition A : B of the sites of 
lattice £ into two subsets A and B, the Schmidt decom- 
position of state 1^1') according to this bipartition reads 

X{A:B) 
a=l 

where Pa, I^'q) and l^*^) appear in the spectral decom- 
position of the reduced density matrices (cf. Eq. l6| 



Al pB 



and where the number of terms x(A : B) in the decom- 
position, known as the Schmidt rank, can be used as a 
measure of entanglement between blocks A and B 

In Ref. [2T] a canonical form for the TTN was pro- 
posed, where each bond index of the TTN corresponds 
to a Schmidt decomposition. That is, in its canonical 
form, a TTN can be regarded as a collection of Schmidt 
decompositions of a state according to a family of nested 
bipartitions A: B oi the system. 

In this work we do not use the canonical form of a 
TTN. However, the use of isometric tensors implies that 
the rank of each bond index in our TTN is given by the 
Schmidt rank x(^ : B) of the corresponding partition. In 



TTN for 10 X 10 Lattice 




particular, the bond dimension x ™ Eq. 13 corresponds 



FIG. 5: Isometric TTN for lattices of 6 x 6, 8 x 8 andlO x 10 
sites as used in the manuscript for the purpose of benchmark- 
ing the performance of the algorithm of Sect. |IV| Notice that 
all TTN have the same structure on the two top layers of 
isometries, whose manipulation dominates the computational 
cost of the algorithm, while they differ in the lower layers. In 
particular, in the 10 x 10 lattice two lower layers of isome- 
tries are required, since a single layer of isometries mapping 
a block of 5 X 5 sites directly into a single effective site would 
have been too expensive given the present capabilities of the 
desktop computers used for the simulations. 



to the Schmidt rank between two halves of the system, 
as well as between one fourth and three fourths of the 
system. 



G. Symmetries 

The symmetries of a state |^) of the lattice C can 
often be incorporated to some extent into the TTN, re- 
sulting in a reduction on computational costs. One can 
distinguish between space symmetries, such as invariance 
under translations e.g. by one lattice site or invariance 
under rotation of the lattice by e.g. 90°, and internal 
symmetries, such as particle number conservation or spin 
isotropy. 

The coarse-graining implicit in the TTN ansatz is in- 
compatible with most space symmetries. As a result, 
a TTN approximation to a symmetric state typically 
breaks such symmetries. However, the symmetry is seen 
to be restored in the limit of a large x- In addition, 
the isometries can often explicitly incorporate part of the 
symmetry. For instance, in approximating states that are 
invariant under translations in 4 x 4, 6 x 6 or 8 x 8 lat- 
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tices by using the TTNs of Figs. [T]and[5j one can choose 
all the isometrics on a given layer of the TTN to be the 
same. 

In contrast, internal symmetries can be implemented 
exactly in the TTN. Suppose for example that the state 
is known to have a well defined particle number (U(l) 
symmetry) or to be a singlet under spin rotations (SU(2) 
symmetry). Then one can choose all the isometrics of 
the tree to be covariant under the action of the symme- 
try, in such a way that: («) the symmetry is preserved 
exactly by any value of x, and (ii) many parameters of 
the isometrics are fixed by the symmetry, leading to a 
significant reduction in computational cost. We refer to 
[60t \6T\ [62] for more details. In the actual computations 
presented in this work we have not implemented internal 
symmetries in the TTN. 



H. Relation to real-space RG 

Being based on coarse-graining the lattice C, the 
present approach is closely related to the real-space RG 
ideas and methods proposed by Kadanofi^, Migdal and 
Wilson H HSl mi UHl [53]. The TTN ansatz can 
indeed be regarded as a specific implementation of the 
spin-blocking schemes that these authors put forward. 

However, it is important to emphasise the differences 
between the present approach and those usually associ- 
ated to real-space RG methods. First of all, here we at- 
tempt to obtain a quasi-exact description of the ground 
state |\E'gs) of a finite lattice C, which forces us to con- 
sider effective sites with a dimension Xt that grows (ex- 
ponentially!) with the number of iterations r of the 
coarse-graining transformation. Instead, real-space RG 
approaches typically attempt to identify and characterise 
the fixed points of the RG fiow on an infinite system and 
consider a fixed dimension Xt- A second important dif- 
ference is in the way the isometrics are chosen. Wilson 
proposed to consider the restriction Hb of the Hamilto- 
nian ff on a block of sites B, and to choose the isome- 
trics in order to preserve the subspace corresponding to 
the lowest energy eigenvalues of Hb- Here, instead, we 
aim at globally minimising H (see Sect. 



IV I, thereby 



following the path initiated with White's density matrix 
renormahsation group (DMRG) [SSI 1571 [55] . 



Relation to MERA 



(i) 



(iv) 





(V) 

o 



FIG. 6: Computation of the expectation value {'I'jo'''' |'I') of a 
one-site operator o'°' acting on site s £ C. (i) Tensor network 
to be contracted, (ii) Tensor network left after many of the 
isometries are annihilated by their Hermitian conjugate, see 
Fig. [2j After the steps from (iii) to (v) the expectation value 
is obtained. 




FIG. 7: Computation of the expectation value (*|o['''o[''''|*) 
corresponding to a two-site correlation function, (i) Tensor 
network to be contracted, (ii) Tensor network left after sev- 
eral isometries are annihilated by their Hermitian conjugate. 
After the steps from (iii) to (v) the expectation value is ob- 
tained. 



The TTN ansatz considered in this work can also be 
thought of as a particular case of the MERA, see Refs. 
[HI [131 [Ml [E] . Specifically, a MERA where disentanglers 
are replaced with identity operators becomes a TTN. We 
emphasise that the manipulations of a TTN, as discussed 
in Sects. Ill and IV, are different than those for the 
MERA. Indeed, the absence of disentanglers changes the 
optimal pattern of tensor network contractions. As a 
result, for instance, the scaling of the computational cost 



with X is much smaller with a TTN (namely C(x'*)) than 
with a MERA [T3]. Of course, for large 2D lattices the 
benefits of having a cost that scales as a smaller power 
of X are offset by the need to use a much larger value of 
X, and simulations with a TTN become more expensive 
than with the MERA. 
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one-site operators Oq' and o\ 



FIG. 8: Computation of the overlap or fidelity {^'i|4'2) be- 
tween two states and 1^2) each represented with a TTN. 
(i) Tensor network corresponding to (*I'i|^2)- Notice that in 
this case no isometry is annihilated, since the isometries of 
the two TTNs are not the same. After the steps from (ii) to 
(iv) the the overlap or fidelity is obtained. 



III. COMPUTATION OF LOCAL OPERATORS, 
FIDELITIES AND ENTROPIES 

In this section we assume that a TTN for the state \^) 
of an Lx L lattice C has been provided, and explain how 
to extract a number of quantities of interest from it. 

This section is presented before explaining the opti- 
misation algorithm in the next section mostly for two 
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reasons. On the one hand, the algorithm of Sect, 
is only one of many possible ways of obtaining a TTN 
(one could alternatively consider using a different opti- 
misation algorithm [H HH [Ml HD , or obtain a TTN 
through an analytical derivation [THJ \T7\ ) and yet in all 
cases it is still necessary to extract information from the 
TTN representation. On the other hand, by explaining 
now how to compute quantities of interest from a TTN, 
we also introduce material that will be useful later on in 
order to understand the optimisation algorithm. 



Expectation value of local operators, two-point 
correlation functions and fidelity 



(17) 



Therefore, without loss of generality we can concentrate 
on the calculation of a two-point correlation function 



mo 



(18) 



This computation is achieved by contracting the tensor 
network of Fig. [Tj A minor difference with the previous 
contraction for a single-site operator is that now less pairs 
of isometries are annihilated. 

More generally, the expectation value of a product of 
p one-site operators (^'lo^^^lol"^! • • • ol^pl j^E') can also be 
obtained by contracting a similar tensor network, and 
so can the overlap or fidelity (^i|^2) between two states 
l^*!) and |^'2) represented by a TTN with equivalent tree 
structure, see Fig. |8] 



B. Spectrum and entropy of the reduced density 
matrix of a block of sites 



Finally, from the TTN it is straightforward to compute 
the spectrum {pa\ of the reduced density matrix p"^ (cf. 
Eq. [6]) for certain blocks A of sites, namely those that 
correspond to an effective site of any of the coarse-grained 
lattices £1, • • • , Ct-i- Fig. [9] illustrates the tensor net- 
work corresponding to p"^ for the case when A is one half 
of the lattice. As before, many pairs of isometries are 
annihilated. In addition, the isometries contained within 
region A can be removed since they do not affect the 
spectrum of From the spectrum {pa}, we can now 
obtain the entropy S{A) of Eq. [t] 

In the large x regime, where the bond dimension at the 
top layers of the TTN is much larger than in the lowest 
layers, the cost of contracting any of the tensor networks 
in Fig. |6][9]is dominated by matrix multiplications whose 
computational cost scales as x^. Thus, this is the cost of 
all the tasks discussed in this section. 



We start by noticing that since the TTN is made of 
isometries, the state \^) it represents is automatically 
normalised, (^l^*) — 1. 

Given a local operator o^^^ that acts on a single site s 
of £, the expectation value 

(*|oW|*) (16) 

can be computed by contracting the tensor network of 
Fig. |6] Notice that an important fraction of the isome- 
tries in the TTN are annihilated by their Hermitian con- 
jugate pair, see Fig. [2] and are therefore not involved in 
the computation of (^'jol^l l^*). 

A local operator o'*" 1 that acts on two sites s and s' 
of C can always be decomposed as a sum of products of 



IV. OPTIMISATION ALGORITHM 

In this section we describe an algorithm to optimise 
the TTN ansatz so that it approximates the ground state 
I^'gs) of a Hamiltonian iJ, 

H\^as)=E^,\^^^), (19) 

defined on an i x i lattice C with torus topology. For 
simplicity we will assume that the Hamiltonian H de- 
composes into two-site terms that couple only pairs of 
nearest neighbour sites s, s' € C, 

H^Y^ /iI"'"'!, (20) 

{s,s') 



htnth ritiith rrtiith 





(iii) 

r 



FIG. 9: Computation of the spectrum {pa} of the reduced 
density matrix for a block A that corresponds to one of 
the coarse-grained sites, (i) Tensor network corresponding to 
p'* where A is half of the lattice, (ii) Tensor network left 
after several isometries are annihilated with their Ifermitian 
conjugate, (iii) since the spectrum of p''*' is not changed 
by the isometries acting on A, we can also eliminate those 
isometries and we are left with a network consisting of only 
two tensors, which can now be contracted together. 



although much more complicated Hamiltonians (e.g. 
with plaquette interactions or arbitrarily long-range in- 
teractions) can be also considered with only minor mod- 
ifications of the algorithm. 



A. Cost function and optimisation strategy 




FIG. 10: Tensor network corresponding to the cost function 
E{{wi}) — {^i^.y\H\'^iiu.y) to be minimised. 




FIG. 11: Tensor network representation for the cost function 
E{w) = F{w) -f G in Eq. 25 depending only on one isometry 



Given the TTN ansatz at a fixed value of Xi 

our goal is to minimise the expectation value 



(21) 



as represented in Fig. [TO] by optimising all the isometries 
{wi} in the TTN, so as to obtain an optimal approxima- 
tion E{{wi}) to the ground state energy Eos, 



min {^-i^^jliJl^-i^^}) 



(22) 



as well as an optimal TTN approximation to the 

ground state I^'gs)- 

An exact solution to Eq. [22] is not known. However, 
one may attempt to approximately minimise the energy 
E{{wi}) in many different ways. Here we will do so by 
means of an iterative optimisation strategy, which is an 
adaptation to the present context of the algorithm de- 
scribed in Ref. [T5]. 

Starting with some set of initially random isometries 
{wi, W2,W3, ■ ■ ■}, we will first optimise one of them, say 
Wi, to obtain an optimal w'l. Then, given the updated set 
{w'i,'W2-,W3, ■ ■ ■}, we will optimise another isometry, say 
W2, obtaining Wj. In the next step, given the updated set 
{w[, 1^3, • • ■ }, yet another isometry will be optimised, 
and so on, until we have optimised all the isometries in 
the TTN. This defines one sweep. Then the sweep is iter- 
ated a number of times, until the cost function E{{wi}) 
is seen to converge according to some criterion, for in- 
stance until it does not change between sweeps by more 
than some small amount. 



B. Optimisation of an isometry 

Next we explain how, given a set of isometries {wi} 
for the TTN at some stage of the minimisation proce- 
dure, one can optimise one isometry w. Recall that w is 
associated to a block A of sites of £. 

First we notice that the cost function E{{wi}) decom- 
poses as a sum of two-site contributions 

m^^.}) = E (23) 

{s,s') 

^ E (24) 

From now on, we also assume for simplicity in the expla- 
nation that h^'^''^ 1 is the product of two one-site operators. 
[If it is not, we can always decompose ft.!^'* ' as a sum of 
such products.] 

When viewed as a function of w only. Fig. [TT] E{{wi}) 
can be divided into two pieces. 



E{w) =F{w) + G. 



(25) 



Here F{w) collects all two-site contributions E^'^ where 
at least one of the two sites s,s & £' are included in the 
block A associated to w, and F collects the rest of two- 
site contributions, see Fig. 12 Notice that if both s and 
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FIG. 12: Examples of the three different types of two-site 
terms E"" contributing to E{w): in (i) both s and s' are 
contained within the block A associated to w; in (ii) only one 
of the sites, say s, belongs to A; finally in (iii) both sites s and 
' are outside A. The terms (i) and (ii) contribute to 



and F^^ (w) in Eq. 26 respectively, whereas the term (iii) 



contributes to the constant G in 1251 





FIG. 13: Examples of the two different types of two-site terms 
that contribute to the environment T for the isometry w: (i) 
both s and s' are contained within the block A associated 
to w, and therefore this term contributes to T'*'*; (ii) only 
one of the sites, say s, belongs to A and therefore this term 
contributes to T'^^. 



s' in lie outside the block A, then the pair w and 



cancels out due to Eq. [2] and E'"' does not depend 
on w. Therefore G is independent of w and we can focus 
on minimising F(w). In turn, F(w) can also be divided 
into two pieces. 



where F^^{w) collects all contributions E'^ 



(26) 
with both 



sites s and s' in A, whereas F^^{w) corresponds to the 
terms E^^ that include one site in A and the other site 
in its complementary B (cases (i) and (ii) of Fig 
The optimisation 



121 



minF(w) 



(27) 



is bilinear in w, and is subject to the isometric con- 
straint of Eq. [2j Unfortunately, once more we do not 
know how to solve this minimisation exactly. 

Following Ref. [15) . we will approximately minimise 
F(w) as follows. First we linearise F{w) by considering 
w to be independent of w\ and then we minimise the 
resulting cost function I{w) — tr(Tu'), 



\.I{vj) — mintr(Tw), 



(28) 



where T is the environment of w. The function I{w) 
can be minimised exactly, with the optimal solution cor- 
responding to w' = —VU\ where T = USV^ is the 
singular value decomposition of T. 

Once we have obtained the optimal w' , we can replace 
w'^ with w''^ in F{w), resulting in an updated environ- 
ment T' that we use to minimise I{w) again. Iteration 
produces a sequence of isometries {w, w', w" , ■ ■ ■} that 
typically lead to monotonically decreasing values of the 
cost function, that is F{w) > F{w') > F{w") > • • • . One 
could in principle iterate the minimisation of F^w) until 
some level of convergence has been reached. However, in 
practice we only use a small number of iterations (even 
just one) before moving to optimise another isometry of 
the TTN, since in actual simulations this is seen to be 
already enough to perform the minimisation of Eq. |22[ 
The order in which the isometries are optimised does not 
seem to play a relevant role in the performance of the 
algorithm. 

All that is left is to explain how to compute the en- 
vironment T of an isometry. Again, the environment 
breaks into two-site contributions corresponding to the 
terms E^" that appear in F'^^l'w) and F"^^{w), 



r^AA 



-AB 



(29) 



Fig 



13 



shows examples of two-site contributions to T'^^ 
and f^^. 

The cost of optimising an isometry comes from the 
computation of the environment T and from its singular 
value decomposition. These costs depend on which isom- 
etry w is optimised, but the cost of sweeping over all the 
isometries of a given layer of the TTN can be seen to 
scale as 0{Lx'^), since there are 0{L) Hamiltonian terms 
/iI"* 1 at the boundary between two halves of the system 
and computing the associated contribution to an envi- 
ronment T has a cost ■ [Notice that the singular value 
decomposition of the environments T for the two largest 
isometries also costs x*]- Therefore the leading order (in 
x) of the cost of sweeping over the whole tree scales as 
0{Lx^), and is also proportional to the number of layers 



10 



in the TTN. In a translation invariant setting where all 
the isometries in a layer of the TTN can be chosen to 
be the same, this leading cost still scales as 0{Lx^) and 
remains proportional to the number of layers, but it has 
a reduced multiplicative pre-factor. In a lattice of 8 x 8 
sites, a typical computation of the ground state (where 
optimisation of isometries proceeds until the expectation 
value of the energy does not change by more than 10^^'^ 
per sweep) requires of the order of 100-500 sweeps. 

In the case of a Hamiltonian made of long range two- 
site interactions, where all sites interact with all sites, 
the number of Hamiltonian terms /i^^^l grows as O(L^) 
and, correspondingly, the overall cost of sweeping over 
the whole tree scales as 0{L'^x^)- 




V. BENCHMARK RESULTS 

In order to test the usefulness of the TTN ansatz and 
to benchmark the performance of the optimisation algo- 
rithm, we consider the quantum Ising model with trans- 
verse magnetic field, as given by the Hamiltonian 

i^Ising = -Y. ® '^ij''' - ^'^''^ ' (30) 

<ss'> 

where Ux and Uz are Pauli matrices and A is the mag- 
nitude of the transverse magnetic field. We consider a 
square lattices made of L x L sites and with toroidal 
boundary conditions. Since each site corresponds to a 
spin-1/2 degree of freedom, its vector space dimension is 
d = 2. In the thermodynamic limit, the model is known 
to undergo a quantum phase transition at a value of the 
transverse magnetic field Ac ~ 3.044 [551 164j . 

We have computed TTN approximations to the 
ground state of i?ising in lattices of linear size L = 
{4, 6, 8, 10, 16, 32} and for several values of x < 500. Ex- 
ploiting translation invariance, we have chosen, on each 
layer of the TTN, all isometries to be the same. For 
i < 8 we are in a quasi-exact regime where results ap- 
pear to be very accurate, whereas for L > 10 we are in 
an approximate regime where the results are not yet con- 
verged with respect to x, but it is still possible to obtain 
qualitatively correct results. 



FIG. 14: Expectation value {a^a^) as a function of the trans- 
verse magnetic field A and for lattices of 4 x 4, 6 x 6 and 8x8 
sites. Notice that, as the lattice size grows, {a^ax) becomes 
steeper and less smooth around A ~ 3, consistent with the ex- 
istence of a critical point at Ac ~ 3.044 in the thermodynamic 
limit [631 164] • 
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FIG. 15: Expectation value (cr^) for the transverse magneti- 
sation as a function of the transverse magnetic field A and for 
lattices of 4 x 4, 6 x 6 and 8x8 sites. Again, as the lattice size 
grows {(Tz) becomes steeper and less smooth around A ~ 3, 
consistent with the existence of a critical point at Ac ~ 3.044 
in the thermodynamic limit [631 164) . 



A. Quasi-exact regime 



For L = 4,6,8 we have computed approximations to 
the ground state of ^^ising for values of the transverse 
magnetic fields in the range A G [1, 5] 



Figs. 



14 



and 



15 



show the sxpectation values for the interaction per link 



1 

2N 



and the transverse magnetisation per site 



(31) 



(32) 



in terms of which the energy per site reads 

-2((Tj.(Ta;) - A((T^) 



(33) 



Both observables in Figs. [T4| and [T5| around A ~ 3 have a 
very steep dependence on A. This behaviour is consistent 
with the presence of a phase transition at A = 3.044 ^ as 



^ In order to find the precise location of the transition point one 
should perform a finite size scaling analysis. However, at least in 
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FIG. 16: Error of the ground state energy per site Ae = e(x) — 
fiei^act (cf. Eq. 33 1 for a 6 x 6 lattice plotted as a function of 
1/x at A = 3.05266. x varies in the range from x = 100 to x = 
550. Ceiact is extracted from Ref. [1]. In the plot we see that 
Ae ranges from Ae ~ 0(10""*) in the case of x = 100 to Ae ~ 
©(lO"'') in the case of x = 550. All the TTN simulations used 
here require much smaller computational resources than the 
full exact diagonalisation calculation as we explain in detail 
in the main text. 



determined in Refs. [631 164]. 

In order to assess the accuracy of our numerical results, 
we first compare the expectation value of the energy with 
its exact value, as obtained in Ref.[T] using exact diag- 
onalisation techniques on a 6 x 6 lattice. [Notice that 
6 X 6 is the largest L x L lattice that can be addressed 
with exact diagonalisation]. For A = 3.05266 (and thus 
near the critical point) the exact value of the energy per 
site as obtained in is e = -3.24727439758.... Fig. 
16 shows the error in the energy per site obtained with a 
TTN with X ranging from 100 to 550. This error is of the 
order of lO^** for x = 100 and under 3x10^^ for x = 550. 
In the first case, the computation lasts 20 minutes on a 
standard PC and uses less than 0.5 Gb of RAM. In the 
latter case, the simulation takes around 2 days and uses 
about 8 Gb of RAM. By comparison, the computation 
by exact diagonalisation required 35 Gb of RAM [IJ. 

This means that, for the model we are considering, we 
obtain accurate results with a fraction of the resources 
needed by the exact diagonalisation algorithms. 

For larger lattices we do not have exact results to com- 
pare against. In this case we study the convergence in 
X of the energy per site e for a value A = 3.05 of trans- 
verse magnetic field close to the critical value Ac. This 



its simplest version, [Tl 1651 . a finite-size scaling analysis requires 
requires knowledge of the energy gap between the ground state 
and first excited state of H. While this gap can in principle be 
computed with a TTN, such computation is beyond the scope of 
the present work, restricted to ground states. 
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FIG. 17: Approximate ground state energy per site e (cf. Eq. 
33 1 for A = 3.05 plotted as a function of 1/x, for lattices of 
6x6 and 8x8 sites. Notice that in both cases the results seem 
to have converged for large x up to several digits of accuracy. 
The insets show the difference Ae = e(x) — e(xmax), as a 
function of 1/x. 



regime is the hardest to simulate, since ground states 
are most entangled at criticality. As shown in Fig. [iTj 
where we plot the energy per site e and its deviation 
Ae = e(x) - e(xmax) from our best estimate e(xmax), 
for values of x around 500, e depends only very weakly 
on X- The figure also shows that, as expected, the 6x6 
case converges faster with large x than the 8x8 case. 
Notice that, bigger systems have higher energies per site, 
in agreement with previous finite size scaling studies [1] 

Further evidence in favour of convergence of the re- 
sults in an 8 X 8 lattice with x is obtained by studying 
the spectrum of the reduced density matrix for one half 
of the lattice. In Fig. 18 we have plotted the largest 200 

3.05. 



eigenvalues of this spectrum, again for A = 3.05. We see 
that in changing x from 200 to 500 in our energy optimi- 
sation, the upper part of the spectrum remains essentially 
unchanged. Also, the spectrum {pa} decays very fast 
as a function of a presenting only around 50 eigenvalues 
larger than 10~^. This also implies that typical errors in 
the expectation value of observables should be very small 
(we will provide examples of this statement in the section 
on the entropies). The study of the spectrum of one half 
of the lattice as a function of A, as displayed in Fig. [19] 
confirms that the ground state is most entangled, and 
therefore its computation most challenging, for A around 
Ac. It is also interesting to notice that, for magnetic fields 
A smaller than the critical Ac, the spectrum presents a 
very peculiar plateaux structure that will be analysed in 
detail in the next section. 

The structure of the TTN manifestly breaks transla- 
tion and rotation invariance and it is natural to ask to 
what degree this affects the structure of correlations in 
the ansatz. Fig. 20 shows the two-point correlation func- 
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FIG. 18: Spectum {pa} for the reduced density matrix p of 
one half of the lattice. The results the ground state of -ffising 
for A = 3.05 in a 8 X 8 lattice. Notice the relatively fast decay 
of the spectrum, with e.g. Pa < 10~* for a > 50. Also, 
calculations with x = 200 and x ~ 500 produce spectra that 
are very similar for small a. This is an indication that the 
largest eigenvalues {pc, for small a) are aheady very close to 
their exact value. 




FIG. 19: Spectrum of the reduced density matrix for one half 
the 8x8 lattice for different values of A. The calculations, con- 
ducted with X ~ 100, show that the spectrum decays slowest 
for A near Ac. It also shows that for magnetic fields smaller 
than Ac, the spectrum develops a clear structure of plateaux. 
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FIG. 20: Two-point correlation functions C{x,y) of Eq. |34] 
for the ground state of iifising with transverse magnetic field 
A — 3.05. The correlation function between distant point in 
the torus remain large, as one would expect of a system that 
becomes critical in the thermodynamic limit. The results, 
obtained with a TTN with x ~ 100, show that the invariance 
of the system under 90° rotations is preserved, in spite of 
the fact that the TTN manifestly breaks it at its top layers. 
Indeed, one can hardly distinguish C(r, 0) from C(0, r). 



The results in this section demonstrate that, for the 
model under consideration, the TTN approach offers a 
reliable route, based on exploiting the entropic area law, 
to extend the domain of quasi-exact results well beyond 
what is possible using exact diagonalisation techniques 



tion 

(34) 

where (x, y) is a vector of integers indicating the position 
of a lattice site. Results obtained for a 8 x 8 lattice with 
just X = 100 hardly show any difference between the cor- 
relation functions in the x and y directions. This seems 
to indicate that the space symmetries expected in the 
ground state have already been restored at a relatively 
small value of x- 



B. Approximate regime 

For lattices of linear size L > 10 we no longer obtain 
convincingly converged results for x ~ 500 when trying 
to approximate the ground state of -ffjsing for A close 
to Ac- Interestingly, however, we still obtain reasonably 
converged results for a large range of A away from Ac, 
which in the case of a 10 x 10 lattice allows us to obtain 
qualitatively the whole phase diagram of the system, see 
Figs. |2l]and|22| 
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FIG. 21: Expectation value {a^Gx) as a function of the trans- 
verse magnetic field A and for a 10 x 10 lattice. The inset 
shows results obtained with x = 100 and x = 200 for values 
of the transverse magnetic field A close to Ac. In this approx- 
imate regime, the TTN algorithm produces results that are 
not converged with respect to x near the quantum critical 
point. 



0.9 

0.8 
0.7 
0.6 




FIG. 22: Expectation value {a^) as a function of the trans- 
verse magnetic field A and for a 10 x 10 lattice. The inset 
shows results obtained with x = 100 and x ~ 200 for values 
of the transverse magnetic field A close to Ac. 



More generally, we find that converged results for lat- 
tices as large as L = 16 and i = 32 can be obtained, with 
X < 500, for values of A not too distant from Ac. This 
can be explained by the presence of a plateaux structure 
in the spectrum of the reduced density matrix of one half 
of the lattice, see Fig. [23] The first plateau consists of 
exactly 2L eigenvalues Pa, that is a e [2,2i + 1]. The 
second plateau is much larger, but its eigenvalues are of- 
ten already very small. For instance, for A = 2.4, the first 
plateau corresponds to Pa ~ 10~^, whereas in the second 
plateau to pa ~ 10~^ — 10" 



FIG. 23: Spectrum {pa} of the reduced density matrix for one 
half of a L X L lattice in the ground state of i/ising for A = 2.4. 
These results, obtained with only x = 100, show the presence 
of a plateau of exactly 2L eigenvalues Pa, separated by two 
or more orders of magnitude from those of the next plateau. 
The structure of plateaux can be understood as a perturbative 
version of the entropic area law and explains why a TTN with 
relatively small x can still produce converged results away 
from the critical point for large lattices L « 10 — 30. For a 
related discussion see also the subsection B of section III in 
Ref. Hg. 




Importantly, Fig. 24 shows 



FIG. 24: Spectrum {pa} of one half of a 10 x 10 lattice in 
the ground state of -ffising for A = 2.4. Results obtained with 
X = 32, 64 and 100 do not differ significantly in the first 21 
eigenvalues. This shows that the presence and composition of 
the first plateau of 2L eigenvalues of Fig. 23 (with L = 10 in 
this case) is robust with respect to x- 



that simulations with a value of x slightly above 2L can 
already accurately reproduce the first plateau and ob- 
tain a reasonable approximation to the ground state of 
the system. This can be explained using perturbation 
theory as done e.g. in Sect. III.B in Ref. [1^ . 



14 



1.4 



1.2 



X 10 



0) 

< 



0.8 
0.6 
0.4 
0.2 













1 \ A e =1.3 10"^- 








Ae=10"^ 


_ 1 






A e =0.5 1 0"- 




4x4 \ 


8x8 


\ 10x10 

^ Ae=10"'*_ 


1* ' ' ~~^^^^T=*===^ , 



100 200 300 400 500 

X 




FIG. 25: Plot of Ae defined in equation |35| for torus of sizes 
L — {4, 6, 8, 10} and various x- We present several choices of 
e in the range 10~* < e < 1-3 10^'^ represented by horizontal 
lines of different colours. The arrows identify the Xe defined 
in the text for the particular choice of £ = 10"'^. 



FIG. 26: Plot of the logarithm of Xe as a function of L for 
several choices of e in the range 10"* < e < 1.3 10"^ In each 
case, for L large enough, the data lie on a straight line con- 
firming that the cost of the simulation increases exponentially 
with the size L of the system. 



C. The exponential cost 

We have seen that the presence of plateaux in the spec- 
tra of reduced density matrices can reduce the cost for 
an approximate description of ground states with a TTN 
to a linear function in the size of the system. In gen- 
eral, however, such short plateaux are not expected (see 
e.g. the valence bond crystal example given in section 



II C). In particular, close to the critical point, we do not 



observe any plateaux. We have repeatedly stated that a 
faithful representation of the ground state with a TTN in 
this regime requires an exponential cost in the size of the 
system. Here we make this statement more precise. In 
order to achieve this we study how the rank of the TTN 
X should increase to keep the error in the energy (as an 
example of a local observable) constant as we increase the 
system size. The error in the energy is estimated from 
our numerical data as 



(35) 



where exma^ our best available result. 

For each system size L we denote by Xe(i) the min- 
imum X that leads to at most an error e in the energy. 
This procedure is illustrated in Fig. 25 In turn. Fig. 



26 displays the value of XeiL) as a function of L. It 
shows that XeiL) grows exponentially with L for large L, 
thereby confirming that the cost of faithfully represent- 
ing the ground state with a TTN increases exponentially 
with the linear size of the system. 



D. Error analysis 

The TTN is a variational ansatz whose precision can 
be improved by increasing the value of the refinement 



parameter x- Ideally, for very large x, the results should 
become exact. A precise theory on how to detect this 
asymptotic regime is beyond the scope of this work. How- 
ever we see that several observables O = (^'^|0|^^) 
converge to their large x value in a characteristic way. 
Namely, 0{l/x) is a monotonic function of 1/x with pos- 
itive, monotonically increasing derivative. For such ob- 
servables, a rough estimate of the error induced by using 
a finite value of x can be obtained as follows. The mono- 
tonic nature of O ensures that 0{\ / Xmax) > C(0) = 
O exact- With a linear fit to the behaviour of 0(l/x) 
close to I / Xmax wc Can extrapolate O to 1/x = and 
obtain Oiow Indeed, the conjectured properties of the 
derivative of O ensure that Oiow ^ C^exact- In this way 
we manage to bound the exact solution with data avail- 
able from the numerical simulations 



Oiow < O exact < C(l/Xmax)- 



(36) 



As important examples we consider in Fig. [27] the be- 
haviour of both the ground state energy and (minus) the 
entanglement entropy of the reduced density matrix of 
half the torus, for the critical Ising Model on a 4 x 4 lat- 
tice. There, we can rely on exact diagonalisation results 
and check that [©(l/xmax) — Oiow\ is an upper bound 
to the error induced by considering smaller x than the 
one required by the exact solution. In this way we have 
estimated the errors appearing in the following section. 



VI. APPLICATION: ENTROPIC AREA LAW 

The study of the entanglement entropy for the ground 
states of 2D quantum systems has been the subject of 
several recent works gZlllHlinilMlEIllSllSllMl- The 
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FIG. 27: Analysis of — S1/2 (top plot) , minus the entangle- 
ment entropy of half a torus, for the critical Ising model on 
a 4 X 4 lattice, as a function of 1/x and the ground state 
energy (bottom plot). The exact result requires x ~ 256. 
If we are only able to compute both observables for x 32 
(the region in the plots on the right of the vertical black line) 
the values of both observables at x = 32 provide an upper 
bound to their exact results. A lower bound is obtained by 
linearly extrapolating the two observables from the two best 
values at X = 16 and x = 32 to 1/x = 0. This is repre- 
sented in the plots by the lower red line on the plots where 
the extrapolated values are denoted by Siow and Eiow In 
this way it is clear that —Siow < —Sexact < — S(1/xt7 
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and Elow < Eexact < E{l/xmax) as expressed in Eq 
in the main text. The differences \S{l/xmax) — Siom\ and 
\E{l/xmax) — Eiow\ thus provlde upper bounds to the error 
for both E and S induced by considering a smaller x than the 
one required by the exact solution. 



TTN approach provides a natural scenario for these stud- 
ies. 



pute the entanglement entropy of one half and one quar- 
ter of the lattice at the critical point A = 3.044. Collect- 
ing results from Refs. gT) SHI SHI |5S1 , at a quantum 
critical point, the entanglement entropy S1/2 of half the 
torus, with total boundary 2L, should scale as 



5i/2(L) = s,2L 



S-l 

2L 



■ iQCPi 



(37) 



where Jqcp (see Ref. [IHl [STJ [Slj ) should be a universal 
constant. We can now numerically confirm the validity 
of Eq. 37 and at the same time extract estimates for 



the coefficients in Eq. 37 including 'Jqcp for the Ising 
universality class on a torus. This is is presented in Fig. 
28 where, by performing a fit to the numerical data with 
Eq. 37 we obtain 



si = 0.06701(107), 
7QCP = 0.030(26), 
s_i = -0.02(14), 



(38) 
(39) 
(40) 



ith 
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rjr = 0.003 . The curve described by Eq 
is, hence, a good description of the scaling form of Si~ 
However the results we obtain are compatible with set 
ting s_i to zero. If we do this, and repeat the fit we 
obtain 



si = 0.06722(18) 
7QCP = 0.0250(21), 



(41) 
(42) 



with 



n df ~ 0-'-'02 (slightly lower than the previous 
case). The values for si and ^qcp £^re also compati- 
ble with the ones of the previous fit. Their accuracy is 
however improved by one order of magnitude. These re- 
sults suggest the absence of the correction term s_i in 
the scaling of the entanglement entropy of half a torus. 
It is interesting to notice that Jqcp, that should be uni- 
versal, is positive as predicted in Ref. |5Tj . 

Compihng results from Refs. [HI HHl |49l US [56] , the 
entropy Si ^4 of one quarter of the torus should scale as 



5i/4(i) = si2L- 



s-i 
2L 



So log 2L + const. (43) 



where the presence of a logarithmic term is induced by 
the corners of the square block [371 UHl IS3 ISS] • Compar- 
ing with Eq. [37] and using that the interacting boundary 
is the same for one half and one quarter of the torus, we 
can extract sq from 



Si/4{L) — Si/2{L) — So log 2L + const. 



(44) 



A. Entanglement entropy 

From the TTN approximation to the ground state of 
^^ising on a L X L torus for L = {4, 6, 8, 10} we can com- 



^ In this section of the work the letter refers to the standard 
name given to the sum of the residual in a least square fit, and 
should not be confuse with the truncation parameter of the TTN 
that is called x through all the rest of work. 
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FIG. 28: Entropy Si/2{L) of one half of the torus as a func- 
tion of the hnear size L, corresponding to the ground state of 
-ffiaing with A = 3.044. The results for L ^ 4,6, 8, 10 confirm 
the linear growth predicted in Eq. |37| where no logarithmic 
correction is expected. The results of our study also seem to 
rule out the presence of the term proportional to s_i. In the 
inset we show the results of the fit. The asymmetric error- 
bars, obtained through the analysis outlined in section |VD| 
are so small that they are hardly visible in the plot. 



This study is presented in Fig. 29 The fit to the numer- 
ical data with Eq. [44| produces an estimate 



So = -0.0381(4) 



(45) 



with a corresponding ^-^ ^ = 0.4. This confirms the 
validity of the scaling form of Eq. |43]as well as a negative 
value for sq in agreement with the theory [571 HHl E3 US] • 

As a side observation, this plot exemplifies once again 
the difference between the quasi-exact and the approxi- 
mate regime. The level of approximation we obtain for 
the L = 10 system implies that those results are almost 
useless in extracting logarithmic corrections to the lead- 
ing scaling laws. 



B. Renyi entropies 

From the TTN approximation to the ground state, we 
can also compute all the Renyi entropies Sn 



Sn 



1 



1 -n 



logtip'^ 



< n < oo 



(46) 



A limiting case of 5„ is given by the single copy entan- 
glement [IS], 



= lim Sn 

n — >oo 



(47) 



is expected to have the same scaling form of Eq. 
37 for the entanglement entropy but with different nu- 
merical coefficients [50l l54l. 



2L 



(48) 



FIG. 29: Difference between entropies for one quarter and one 
half of the torus as a function of the linear size L, correspond- 
ing to the ground state of ifiaing with A = 3.044. The results 
for L — 4,6, 8, 10 allow us to confirm the logarithmic depen- 
dence predicted in Eq. |44[ which is attributed to the presence 
of corners in the boundary of our block for one quarter of the 
lattice. 



The results for Ei of half torus are shown in [30] The 
fit to the numerical data for L = 4,6,8, 10 with Eq. [48| 
produces 



ei 

i'qcp 



0.01724(20), 

0.0499(51), 

-0.161(30), 



(49) 
(50) 
(51) 



with 



n.d.f. 



— 0.0005. This reveals a very good agreement 



between Eq. 48 and the numerical data. It is interesting 
to notice that the coefficient e_i is non-zero and neg- 
ative, in agreement with the theory |53j . In addition, 
as already anticipated by the results contained in Ref. 
[50l [54] , the numerical values of the parameters for scal- 
ing of the single copy entanglement and entanglement en- 
tropy are different. Nevertheless, for the universal term, 
we find that 



i'qcp = 27QCP 



(52) 



to our numerical precision. This is reminiscent of a very 
similar result obtained for one dimensional critical chains, 
where the universal coefficient of the logarithmic scaling 
of the entanglement entropy with the size of the interval 
is two times (instead than one half as in Eq. 52 ) the 



analogous coefficient for the scaling of the single copy 
entanglement [67l [68] . 



VII. DISCUSSION 

In this manuscript we have described a numerical tech- 
nique based on a TTN to compute ground state prop- 
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FIG. 30: Single Copy Entanglement E^^^L) for one half of 
the torus as a function of the linear size L, corresponding to 
the ground state of iifising with A — 3.044. The results for 
L = 4,6, 8, 10 confirm the linear growth predicted in Eq. |48[ 
In this case the term e_i, analogous to the term s_i for the 
entanglement entropy, is negative as predicted by the theory. 
In the inset we present the results of the fit. The asymmetric 
errorbars, obtained through the analysis outlined in section 
|VD| are so small that they are hardly visible on the plot. 



erties of 2D lattice systems. The approach exploits the 
entropic area law and has a cost that scales exponentially 
in the linear size of the lattice. Its goals are necessarily 
more modest than those of scalable tensor network algo- 
rithms such as PEPS and MERA t41|5J[Sl|71IHl[Sl[ini[IIl 

For the model we have considered here, the TTN ap- 
proach offers a simple, effective way of obtaining quasi- 
exact results well beyond what is possible with exact di- 
agonalisation techniques [TJ [5]. We expect that similar 
gains would also occur for other models. We envisage 
that this technique will become a useful tool both to 
study small lattice systems and in investigations based 
on finite size scaling. A highlight of the approach is its 
simplicity, specially when compared to the scalable ten- 
sor network algorithms. In addition, it can be used to 
study block entropies, a task that becomes much less 
straightforward with other methods. 

The TTN algorithm is closely related to the DMRG 



algorithm applied to 2D lattices. It is beyond the scope 
of the present work to conduct the detailed analysis re- 
quired to establish how the performances of the two al- 
gorithms compare. Nevertheless, some preliminary ob- 
servations can be made. Updating the matrix product 
state (MPS) used in DMRG has a cost of 0{x^L'^) per 
sweep, while updating the TTN costs 0{x'^L). This al- 
lows DMRG to consider values of x that are about 10 
times larger with similar computational cost. On the 
other hand the TTN has better connectivity. In a TTN 
all lattice sites are connected through the product of at 
most O(logL) tensors. Instead, when an MPS is used to 
encode the ground state of a 2D lattice, nearest neighbour 
lattice sites are typically connected through the product 
of 0{L) tensors, with a fraction of the sites being con- 
nected through the product of 0{L^) tensors (on a torus). 
As a result, we expect convergence to the ground state 
to be faster using a TTN. In addition, space symmetries 
can be (partially) incorporated in a TTN. 

The TTN is particularly fitted to study entropies and 
their scaling with the size of the system. In this work we 
have reported some numerical results that are compatible 
with the expectation drawn from Refs. [47l HHl [49l [55] 
about the presence, in the scaling form of the entangle- 
ment entropy, of both additive logarithmic and constant 
corrections to the area law. Our results suggest the ab- 
sence of a 1/L correction. After the first draft of our 
paper was presented, a systematic study of all the Renyi 
entropies was reported in Ref. 1541 . This motivated u s to 
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also consider the single copy entanglement E'-^\ Eq 
We have confirmed that the scaling of E''^^ also includes 
a constant additive correction to the area law and a term 
proportional to 1/L, with numerical values different from 
the ones present in the scaling of the entanglement en- 
tropy. This is a hint to the presence of a different set 
of universal constants for each of the Renyi entropies as 
stated in Ref. [SI] . Further aspects of this scenarios can 
be found in Ref. [69] . 
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cial support from the Australian research council (APA, 
FF0668731, DP0878830) is acknowledged. 



[1] C. J. Hamer, Journal of Physics A: Mathematical and 

General 33, 6683 (2000), ISSN 0305-4470, URL http^ 

//www . iop . org/EJ/abstract/0305-4470/33/38/303/ 

[2] M. Henkel, Journal of Physics A: Mathematical and Gen- 
eral, 17, L795-L799 (1984). 
[3] A. Laeuchli and C. Lhuillier, 0901.1065 (2009), URL 

http : //arxiv . org/abs/0901 . 1065 
[4] F. Verstraete and J. I. Cirac, cond-mat/0407066 (2004), 
URL http : //arxiv . org/abs/cond-mat/0407066i 



[5] V. Murg, F. Verstraete, and J. I. Cirac, Physical Review 
A (Atomic, Molecular, and Optical Physics) 75, 033605 
(2007), URL http : //link . aps . org/abstract/PRA/v75/. 
I e03 3605 ' 

[6] J. Jordan, R. Orus, G. Vidal, F. Verstraete, and J. I. 
Cirac, Physical Review Letters 101, 250602 (2008), URL 
http : //link . aps . org/abstract/PRL/vl01/e250602 

[7] V. Murg, F. Verstraete, and J. I. Cirac, Physical Review 
B (Condensed Matter and Materials Physics) 79, 195119 



18 



(2009), URL http : //link, aps . org/abstract/PRB/v79/ 



[9 



[11 



[12: 



[13 



[ 



[14 



[15 



[16: 



[T7 



[is: 



r 

[19 

L 



[20: 



[21 



[22 



[23: 

[24 
[25 



[26 

[27 
[28 
[29 



el95119 

G. Sierra and M. A. Mart in-Delg ado, cond-mat/9811170 
(1998), URL jhttp : // airxiv . org/abs/cond-mat/ 
98111701 

N. Maeshima, Y. Hieida, Y. Akutsu, T. Nishino, 
and K. Okunishi, cond-mat/01 01360 (2001), 
_Phys.Rev.E64:016705,2001, URL |htt:p:7/arxiv.| 
org/abs/cond-mat/0101360; 

Y. Nishio, N. Maeshima, A. Gendiar, and T. Nishino 
(2004), 
04011151 



URL 



Z. Gu, M. Levin, and X. Wen, Physical Review B 
(Condensed Ma tter and Materials Physics) 78, 205116 
(2008), URL h ttp : //link, aps . org/abstract/PRB/v78/| 



^205116 

Vidal, 



G. 



Physical Review Letters 99, 220405 



(2007), URL http: //link, aps . org/abstract/PRL/v99/ 



6220405 

G. Vidal, quant-ph/0610099 (2006), Physical Review 
Letters 101, 110501 (2008), URL .http://arxiv.org/, 



abs/quant-ph/0610099| 



G. Evenbly and G. Vidal, Physical Review Letters 102, 

180406 (2009), URL http://link.aps.org/abstract71 

Q PRL/vl02/el80406 

G. Evenbly and G. Vidal, Physical Review B (Condensed 
Matter and Materials Physics) 79, 144108 (2009), URL 
http : //link . aps . org/abstract/PRB/v79/e l44108, 
M. Fannes, B. Nachtergaele, and R. F. Werner, Journal 

_of Statistical Physics 66, 939 (1992), URL jhttp : //d x . 
doi . org/10 . 1007/BF01055710 

H. Niggemann, A. Klmper, and J. Zittartz, Zeitschrift 
fr Physik B Condensed Matter 1 04, 103 (1997), URL 
|http : //dx . doi . org/10 . 1007/s00 2570050425 
B. Friedman, Journal of Physics: Condensed Matter 9, 

_9021 (1997), ISSN 0953-8984, URL http: //www. io^ 
org/EJ/abstract/0953-8984/9/42/016/, 
M. Lepetit, M. Cousy, and G. M. Pastor, European Phys- 
ical Journal B 13, 421 (2000), URL jhttp : //adsabs . | 
harvard . e du/abs/2000EP JB . . .13. .42lO 

J. 



M. A. Martin-Delgado, J. Rodriguez-Laguna, and 
G. Sierra, Physical Review B 65, 155116 (2002), URL 
http : //link . aps . org/abstract/PRB/v65/el55116 
Y. Shi, L. Duan, and G. Vidal, Physical Review A 
(Atomic, Molec ular, and Optical Physics) 74, 022320 
(2006), URL ht tp : //link, aps . org/abstract/PRA/v74/| 
11022320 

D. Nagaj, E. Farhi, J. Goldstone, P. Shor, and 
L Sylvester, Physical Review B (Condensed Matter and 
Materials Physics) 77, 214431 (2008), URL httpT/T] 
J.ink . aps . org/abstract/PRB/v77/e214431 , ^ 
~R. Otsuka, Physical Review B 53, 14004 (1996), URL 
|http : //link . aps . org/abstract /PRB/v53/pl4004 
L. P. Kadanoff, Physics 2, 263 (1966). 
L. P. Kdanoff and et al. Reviews of Modern Physics 39, 
_395 (1967), URL |http : //link . aps . org/abstract /RMP/| 
"v39/ p395 

K. G. Wilson, Reviews of Modern Physics 47, 773 (1975), 
URL http : //link . aps . org/abstract/RMP/v47/p773 
M. E. Fisher, Reviews of Modern Physics 70, 653 (1998), 
URL http : //link . aps . org/abstract /RMP/v70/p653 

T. W. Burkhardt and J. M. van Lceuwcn, Real-space 

renormaUzation (Springer- Verlag, Berlin, 1982). 

L. P. Kadanoff, Physical Review Letters 34, 1005 (1975), 



[30 
[32 



jhttp : // arxiv . org/abs/cond-mat/ [ 



[33 



[34 



[35 
[36 

r 

[37 
[38 

r, 
[39 



[4o: 

[41 

L 



[42: 



[43 



[44 
[45: 



[46: 



[47: 



[48: 

[49 



[52: 



URL http : //link . aps . org/abstract /PRL/v34/pl005 
L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, Phys- 



^cal Review D 34, 373 (1986), URL |http : //link . apsT 
org/abstract/PRD/v34/p373, 

M. Srednicki, Physical Review Letters 71, 666 (1993), 
URL http : //link . aps . org/abstract /PRL/v71/p666 
M. B. Plenio, J. Eisert, J. Dreissig, and M. Cr amer, 
Physical Review Letters 94, 060503 (2005), URL httJTI 



//link . aps . org/abstract/PRL/v94/e060503 

TTMasanes, 0907.4672 (2009), URL http: //a ^iv. org/| 
_abs/0907.4672, 

J. I. Latorre, E. Rico, and G. Vidal, quant-ph/0304098 
_(2003), Quant. Inf. Comput. 4 (2004) 48-92, URL http : | 
//arxiv . org/abs/quant-ph/0304098 
J. Eisert, M. Cramer, and M. B. Plenio, 0808.3773 
(2008), URL http : // arxiv . or g/abs/0808 . 3773, 
S. R. White and R. M. Noack, Physical Review Letters 



68, 3487 (1992), URL http : //link . aps . org/abstract/ 
PRL/v68/p3487, 

S. R. White, Physical Review Letters 69, 2863 (1992), 
URL http : //link . ap s . org/abstract/PRL/v69/p2863 
R. M. Noack and S. R. White, Physical Review B 



47, 9243 (1993), URL jhttp : //link . aps . org/abstract/ 
PRB/v47/p9243, 

U. SchoUwock, Reviews of Modern Physics 77, 259 
(2005), URL http : //link . aps . org/abstract /RMP/v77/, 
•p259l 



N. Shibata, J. Phys. A: Math. Ge n. 36, R381 (2003), 
URL http : // arxiv . org/abs/cond-mat/0310028 

S. R. White and A. L. Chernyshev, Physical Review 
Jl^etters 99, 127004 (2007), URL http:/ /link.aps.org/ 
^stract/PRL/v99/el27004 

S. R. White and D. J. Scalapino, Physical Review Letters 
_80, 1272 (1998), URL ht tp : //link . aps . org/abstractT] 
"PRL/v80/pl272| 

S. Liang and H. Pang, Physical Review B 49, 9214 



(1994), URL |http : //link. aps . org/abstract/PRB/v49/ 
p9214 



P. Hen elius, Physical Review B 60, 9561 (1999), URL 

|http:// link . aps . org/abstract/PRB/v6 0/p9561. , 

M. S. L. du Croo de Jongh and J. M. J. van Leeuwen, 



_Pliysical Review B 57, 8494 (1998), URL [http: //link. 
aps . org/abstract /PRB/v57/p8494 

T. Xiang, J. Lou, and Z. Su, Physical Review B 64, 
104414 (2001), U RL http : //link . aps . org/abstract/] 
PRB/v64/el04414| 

H. Casini and M. Huerta, Nuclear 
Physics B 764, 183 (2007), ISS N 0550- 
3213, URL http: //www. sciencedirectT] 

com/science/article/B6TVC-4MMWHHD-3/2/ j 
;^bfa856cl91dba6b0db26bdlda7eb3cl 
E. Fradkin and J. E. Moore, cond-mat /0605683 
_(2006), Phys.Rev.Lett.97:050404,2006, URL [http:// 
arxiv . org/abs/ cond-mat/0605683 

B. Hsu, M. Mulhgan, E. Fradkin, and E. Kim, Physical 
Review B (Condensed Matter and Materials Physics) 79, 
115421 ( 2009), URL http : //link . aps . org/abstract/, 



PRB/v79/ell5421 



J. Stephan, S. Furukawa, G. Misguich, a nd V. Pasquier, 
0906.1153 (2009), URL |htt p : //ar:!dT. org/abs/0906T] 

T153 

E. Fradkin, 0906.1569 (2009), URL http: //ar xiv. org/] 
abs/0906.1 569 

R. Yu, H. Saleur, and S. Haas, Physical Review B 



19 



(Condensed Matter and Materials Physics) 77, 140402 
(2008), URL http : //link . aps . org/abstract/PRB/v77/, 
el40402l 



[53] H. Casini and M. Huerta, 0905.2562 (2009), URL |http : 
I //ar xiv.org/abs/0905 . 2562 

[54] M.^7 Metlitski, C. A. Fuertes, and S. Sachdev, Physical 
Review B (Condensed Matter and Materials Physics) 80, 
115122 (2009), URL |http : //link . aps ■ org/abstract/| 
Q PRB/v80/ell5122 

[55] H. Casini and M. Huerta, Journal of Statistical 
Mechanics: Theory and Experiment 2005, P12012 
(2005), ISSN 1742-5468, URL http://www.iop.org/EJ/ 
abstract/1742-5468/2005/12/P12012/ ^ 
[56] T. Nishioka, S. Ryu, and T. Takayanagi, 0905.0932 

(2009), URL http://arxiv.org/abs/0905.0932 
[57] C. H. Bennett, H. J. Bernstein, S. Popescu, and B. Schu- 
macher, Physical Review A 53, 2046 (1996), URL http: 



I / /link . aps . org/abstract /PRA/v53/p2046 
[58] L Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys- 

ical Review Letters 59, 799 (1987), URL |http : //link .| 

I ap s . org/abstract/ PRL/v59/p799 , 
[59] G. Vidal, Physi cal Review Letters 91, 147902 



(2003), URL |http : //link. aps . org/abstract/PRL/v91/ 



_el47902 

[60] S. Singh, R. N. C. Pfeifer, and G. Vidal, 0907.2994 

(2009), URL http://arxiv.org/abs/0907.2994 
[61] L P. McCuUoch and M. Gulacsi, EPL (Europhysics Let- 
ters) 57, 852 (2002), ISSN 0295-5075, URL http: //www., 
I iop ■ org /EJ/abstract/0295-5075/57/6/852/ 



[62] S. Singh, H. Q. Zhou, and G. Vidal, cond- 
m at/070142 7 (200 7), URL http : //arxiv . org/abs/, 
I con d-iiiat/0 701427| 

[63] H. Rieger and N. Kawashima, The European Physi- 
cal Journal B - Condensed Matter and Complex Sys- 
tems 9, 233 (1999), URL |http : //dx . doi . org/10 . 1007/| 

I S100510050761 

[64] H. W. J. Blote and Y. Deng, Physical Review E 66, 

066110 (2002), URL http : //link . aps . org/abstract/, 

I PR E/v66/e066110 

[65] M. P. Nightingale, Physica A: Statistical and 
Theoretical Physics 83, 561 (1975), ISSN 
0378-4371, URL http: //www, sciencedirect. 



com/science/article/B6TVG-46CCBJV-85/2/ 



[ 

I bd e4044509a99bc0ee486c9d663048c6l 
[66] J. Eisert and M. Cramer, Physical Review A 72, 042112 
(2005), URL http : //link . aps . org/abstract/PRA/v727] 
e042112 

[67] R. Orus, J. I. Latorre, J. Eisert, and M. Cramer, Physical 
Review A (Atomic, Molecular, and Optical Physics) 73, 
060303 (2006), URL ,http : //link . aps . org/abstract/, 

I PRA/v7 3/e060303l 

[68] I. Peschel and J. Zhao, Journal of Statistical Me- 
chanics: Theory and Experiment 2005, P11002 
(2005), ISSN 1742-5468, URL http: //ww w. iop. org/EJ/| 
abstract/1742-5468/2005/ll/P11002/ 

[69] F. Gliozzi and L. Ta gliacozzo, 0910.3003 (2009), URL 
|http : //arxiv . org/abs/0910 . 3003| 



